%parameter: realizations 
realizations = 25000;

%parameters Xa, Xb: X min and max values
Xa = -0.5;
Xb = 0.5;

%parameters Ya, Yb: Y min and max values
Ya = -0.5;
Yb = 0.5;

%constant: Int. Represents the analytical exact value of the calculation
Int = 0.148496

YBorderTerm = 1 / (Yb - Ya)
XBorderTerm = 1 / (Xb - Xa)

ihat = 0;

for i = 1:1:realizations
	
	x = rand() - 0.5;
	y = rand() - 0.5;
	
	squareSum = (x^2 + y^2);

	f = squareSum * e ^ (- squareSum / 2) / (YBorderTerm * XBorderTerm);
	ihat = ihat + f;
	
	value(i) = ihat / i;
	
	errors(i) = abs(Int - ihat / i);
	
end
media = 0;
for i = 1:1:realizations
	media = media +  value(i);
end

media = media / realizations;
suma = 0;
j = 2;
tenpercent = 0;
while j <= realizations
	valor = (value(j) - media)^2;
	suma = suma + valor;

	desv(j) = sqrt(suma / (j - 1));
	if desv(j)  < 0.1 * Int && tenpercent == 0
		tenpercent = j
	end
	
	j = j + 1;
end

desvio = desv(j-1)
valorfinal = value(j - 1) 

r = 1:realizations;

figure;
plot(r,value);
%title('Integral');
xlabel('Numero de Realizaciones');
ylabel('Integral');


figure;
plot(r,errors);
%title('Error');
xlabel('Número de Realizaciones');
ylabel('Error');

q = 2:realizations;
figure;
plot(q,desv);
%title('desv');
xlabel('Numero de Realizaciones');
ylabel('Desvio');